rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
#8.1 Poisson回归模型
#install.packages("robust")
library(robust)
data(breslow.dat)
str(breslow.dat)
breslow <- breslow.dat[,6:9]

summary(breslow)

hist(breslow$Ysum,
     breaks = 20,
     xlab = "Seizure Count",
     main = "Distribution of Seizure")

#可知并非对称分布
fit <- glm(Ysum~Base+Age+Trt, data = breslow.dat, family = poisson)
summary(fit)
#三个变量皆有统计学意义
#提取模型的系数
coef(fit)
#指数化
exp(coef(fit))
exp(confint(fit))
#汇总
library(epiDisplay)
idr.display(fit)
write.csv(idr.display(fit)$table, file = "poisson.csv")

#8.2 过度离散的判定及处理
#Poisson分布要求均值与方差相等；
#实际过程中，往往出现方差大于均值的过度离散的现象
#遗漏某个重要的预测变量、观测之间彼此不独立

#install.packages("qcc")
library(qcc)
qcc.overdispersion.test(breslow$Ysum, type = "poisson")
#表明存在过度离散
library(epiDisplay)
poisgof(fit)
#拟合优度检验的p值很小，表明模型拟合较差，
#Poisson回归的假设条件不满足，很可能存在过度离散的现象

#1 使用拟Poisson回归
#只需要将glm改为quasipoisson
fit.od <- glm(Ysum~ Base + Age + Trt,
              family = quasipoisson,
              data = breslow)
summary(fit.od)
#回归系数与使用Poisson分布所得相同；但其标准误更大

#2 指定误差项服从负二项分布，使用负二项回归
fit.nb <- glm.nb(Ysum ~ Base + Age + Trt, data = breslow)
summary(fit.nb)
#在考虑过度离散后，变量Age与Trt都没有统计学意义了
poisgof(fit.nb)
#通过了拟合优度检验
AIC(fit,fit.nb)

#8.4 对数线性模型
#Poisson分布用于列联表资料的分析
#不区分因变量与自变量，分析理论频数与实际频数的差异
#对数线性模型为层次模型，如果模型中包含了某几个变量的高阶相互作用项
#这几个变量的低阶相互作用项与主效应项也必须包含在模型中
#一般以饱和模型开始，通过后退法逐步排除没有统计学意义的项，
#最后得到拟合最优的简化模型
#饱和模型为包含了所有变量的主效应、所有低阶相互作用和高阶相互作用的项

#导入数据
dat <- array(c(100,63,4,2,36,84,10,25),
             dim = c(2,2,2),
             dimnames = list(pill = c("no","yes"),
                             gene = c("V-", "V+"),
                             group = c("control","case")))
dat
dat <- as.table(dat)
mydata <- as.data.frame(dat)
mydata
str(mydata)

#建立饱和模型
mod1 <- glm(Freq ~ pill*gene*group, family = poisson, data = mydata)
summary(mod1)
#包含三个主效应、三个一阶交互作用与一个二阶交互作用以及一个常数项

#用step简化模型，标准为AIC值
mod2 <- step(mod1)

summary(mod2)

exp(coef(mod2))

#拟合优度检验
library(epiDisplay)
poisgof(mod2)
#最后表明，模型拟合的效果很好
